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Properties of disordered nematic elastomers and gels are theoretically investigated with empha- 
sis on the roles of non-local elastic interactions and crosshnking conditions. Networks originally 
crosslinked in the isotropic phase lose their long-range orientational order by the action of quenched 
random stresses, which we incorporate into the affine-deformation model of nematic rubber elasticity. 
We present a detailed picture of mechanical quasi- Goldstone modes, which accounts for an almost 
completely soft polydomain-monodomain (P-M) transition under strain as well as a "four-leaf clover" 
pattern in depolarized light scattering intensity. Dynamical relaxation of the domain structure is 
studied using a simple model. The peak wavenumber of the structure factor obeys a power-law-type 
slow kinetics and goes to zero in true mechanical equilibrium. The effect of quenched disorder on 
director fluctuation in the monodomain state is analyzed. The random frozen contribution to the 
fluctuation amplitude dominates the thermal one, at long wavelengths and near the P-M transition 
threshold. We also study networks obtained by crosslinking polydomain nematic polymer melts. 
The memory of initial director configuration acts as correlated and strong quenched disorder, which 
renders the P-M transition non-soft. The spatial distribution of the elastic free energy is strongly 
dehomogenized by external strain, in contrast to the case of isotropically crosslinked networks. 

PACS numbers: 61.30.Cz, 61.41. -He, 64.70.Md 



I. INTRODUCTION 

Elastomers and gels are intrinsically disordered solids 
that retain the memory of their initial states. The non- 
equilibrium nature of their fabrication processes causes 
frozen heterogeneities in the network structure, which 
range in size from mesoscopic to macroscopic scales . 
The presence of the quenched disorder comes to the fore 
when we introduce some soft order in the system. For in- 
stance, density fluctuations of swollen gels near the crit- 
ical point are strongly enhanced by the heterogeneities. 
Under stretching, they produce the so-called "abnormal 
butterfly" pattern in small-angle neutron scattering in- 
tensity It illustrates how the elasticity of gels 
gives rise to a non-trivial effect unexpected in other ran- 
dom systems. Here we address another example of soft 
order in disordered elastic networks. 

Nematic liquid-crystalline elastomers and gels consti- 
tute a unique class of solids characterized by a coupling 
between the orientational and translational degrees of 
freedom. Physical consequences of the strain-orientation 
coupling have been the subject of a considerable amount 
of studies, both theoretical and experimental. Notable 
theoretical advances in the past include: (i) De Gennes |^] 
showed that a spontaneous elongation along the direc- 
tor is induced by the isotropic-nematic (I-N) transition; 
(ii) A molecular model of nematic networks was con- 
structed by Warner et al. extending the classical 
affine-deformation model of rubber elasticity; (iii) Uni- 
formly oriented networks possess soft modes of strain 



and orientation fluctuations that do not accompany any 
change of rubber-elastic free energy. It was first predicted 
by Golubovic and Lubensky ||^ on a phenomenological 
ground and later extended by the afhne-deformation the- 
ory |l^,|ll[. Thus and in other ways, the behavior of ho- 
mogeneous and clean nematic networks is now fairly well 
understood |]l2| . 

In practice, however, nematic networks in equilib- 
rium quite often exhibit polydomain director textures, 
where the orientational correlation length is typically 
in the micron range. Under external strain, polydo- 
main networks undergo a structural change into a macro- 
scopically aligned monodomain state, where the direc- 
tor lies along the extensional direction. This change, 
called the polydomain-monodomain (P-M) transition, 
is characterized by a highly non-linear mechanical re- 

The strain-stress curve 



sponse pj,y,pj,n6 17dJJl 



shows a small slope in the partially aligned (polydomain) 
state. Depending on the material and the method of syn- 
thesis, the slope is sometimes vanishingly small while it is 
sizable in other cases. The macroscopic stress as a func- 
tion of strain shows a steep rise as the system turns into 
the monodomain state. 

There have been a few theoretical attempts to de- 
scribe polydomain networks and their mechanical re- 
sponses. Ten Bosch and Varichon set up the first 
model, in which they attributed the origin of the equi- 
librium texture to a random anchoring field exerted by 
network crosslinks. An interesting analogy with random 
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anisotropy magnets [ pT[ was pursued by Fridrikh and Ter- 
entjev [p2| . They proposed a mapping to the XY model 
under random and homogeneous magnetic fields, from 
analysis of which they predicted a discontinuous stress- 
orientation curve. 

Nonetheless, the role of strain-orientation coupling in 
polydomain networks is still far from clear. There are 
two aspects to be considered. Firstly, the previous theo- 
ries assume only local interactions between domains, for 
instance by arguing that the elastic free energy local- 
izes in domain walls under strain [ p2| . In general, how- 
ever, inhomogeneities in an elastic material cause non- 
local or long-range interactions mediated by the strain 
field. Such elastic interactions control the physics of var 
ious systems, such as solids with dislocations [P3| or sur 
face defects |^,^, phase separating alloys 
undergoing swelling 



7|, gels 

and membranes with inclu- 
sions | |30| |. Disordered nematic networks provide another 
intriguing example, and differ from any of the above ma- 
terials in having a non-scalar order parameter. Secondly, 
the mechanical response should strongly depend on the 
crosslinking condition. Polydomain elastomers have been 
obtained by either of the following ways (i) to 

crosslink a polymer melt in the isotropic phase and then 
cool it into the nematic phase; (ii) to crosslink a nematic 
polymer melt containing polydomain textures. These two 
cases have not been theoretically well distinguished so 
far. We shall refer to them as the cases of isotropic and 
anisotropic crosslinkings, respectively. 

Recently, we studied the elastic interaction in isotrop- 
ically crosslinked networks and found an almost 

completely soft P-M transition The macroscopic 

stress due to the strain-orientation coupling was shown 
to be slightly negative and of 0{a^) in the P-M transition 
regime, where a is the degree of chain anisotropy. This 
contrasts with the earlier prediction of a positive stress of 
0{a) 1^^. The elastic interaction also produces a "four- 
leaf clover" pattern in the depolarized light scattering 
intensity, which resembles the experimental observation 
by Clarke et al. 

In this article, we extend previous work |^,^ and 
provide the details of our picture of the P-M transition. 
Here let us summarize the ideas and results which we 
have not emphasized in previous work. Firstly, we pur- 
sue the idea that random internal stresses destroy the 
long-range orientational order, which was suggested (but 
not proven) earlier in a broader context ||]. This will 
be done by incorporating the notion of frozen heteroge- 
neous strains into the extended affine-deformation the- 
ory We argue that the random internal stresses act 
as stronger sources of disorder than the random molecu- 
lar field due to crossfinks (2^,^. Secondly, evolution of 
domain structure with and without external stretching is 
numerically simulated by a simplified dynamical model. 
The "four-leaf clover" scattering pattern has four peaks 
at finite wavenumbers, and the peak height is a non- 
monotonic function of the macroscopic strain, in quali- 
tative agreement with experiment. We find a slow dy- 



namical relaxation of the structure factor, and show that 
the peak wavenumber asymptotically goes to zero in the 
long-time limit. Thirdly, we study the case of anisotropic 
crosslinking. In this case, the initial director configura- 
tion of a macroscopic polydomain texture is memorized 
into the network. It provides a source of strong and 
correlated disorder, resulting in a non-soft P-M transi- 
tion. The spatial distribution of elastic free energy in 
anisotropically crosslinked networks is strongly dehomog- 
enized by strain, while that in isotropically crosslinked 
networks is unchanged during the P-M transition. 

This paper is organized as follows. In Section II, we 
introduce a random stress model, derive an effective free 
energy, and discuss the mechanism of the soft mechanical 
response. Section III describes a numerical simulation of 
the polydomain state and the P-M transition. In Section 
IV, we analyze the effect of random stresses on director 
fiuctuations in the monodomain state. We study net- 
works prepared in the nematic phase in Section V. In 
Section VI, we summarize the results in comparison to 
existing experiments, and conclude with a proposal of 
future directions. 



II. MODEL AND ANALYSIS 



A. Random stresses in isotropic networks 



It is known since long ago 1^,^ that the network 
structure of gels are often heterogeneous on many length- 
scales, which are considerably larger than the mesh size 
(see Fig.|l|). In the swollen state, these imperfections 
manifest themselves as density inhomogeneity and are 
observed through the so-called butterfiy pattern in neu- 
tron scattering intensity or as speckles in light scattering 
experiment (^J^. Although less discussed in the litera- 
ture, it is natural to expect that elastomers, often fabri- 
cated by drying gels, also contain the memory of hetero- 
geneous network formation. The frozen heterogeneities 
refiect the non-equilibrium nature of the crosslinking pro- 
cesses, and produce random internal stresses in the ma- 
terial. While the roles of random stresses in gels and 
other amorphous solids have been discussed from a phe- 
nomcnological point of view much remains to be 

done to understand them on the basis of a molecular the- 
ory In this subsection, we recapitulate the notion 
of random stresses using the classical affine-deformation 
theory of isotropic rubber networks [^6| , in order to pre- 
pare for modeling disordered nematic networks in the 
next subsection. 

The basic object in the affine-deformation theory is the 
probability distribution of the chain's end-to-end vector 
p. The distribution function at thermal equilibrium is 
isotropic and Gaussian, and given by 



Peq{p) = A/" ^exp 



2f7' 



(1) 
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where = (p^)^^ is a constant, d is the spatial dimen- 
sion, and Af = J dpPeq{p) is the normahzation factor. 
The macroscopic deformation of the network is described 
by the Cauchy deformation tensor. 



A, 



drj 



(2) 



where r and ro are the positions of material points at 
observation and at the moment of crosslinking, respec- 
tively. The basic assumption of the theory is that each 
chain's end-to-end vector affinely changes as p — > A • p 
in response to the macroscopic deformation. The free 
energy per chain is given by pq] 



chain 



,T J dpPnip) InPeq(A-p), 



(3) 



where Po{p) is the probability distribution function at 
the moment of crosslinking, which is not necessarily iden- 
tical to the equilibrium distribution. Wc assume that the 
chains are distorted before crosslinking, and denote the 
deviation from the equilibrium conformation by a tensor 
R, defined by 



{pp),^n{\ + R), 



(4) 



where I is the unit tensor. If the deviation is not very 
large and the chains are not stretched out, we may still 
approximate Pq to be Gaussian, and put 



Po{p) = A/'o ^exp 



2n 



p.(l + R)-i.p 



Substituting (Q) and (||) into Eq.(|3|), we have 



/< 



chain 



TtG + R : G-|-lndet(l + R) 



(5) 



(6) 



where = AkiAkj is the metric tensor of deformation. 
The equation is not new and essentially contained 
in the classical theory of Flory |Q . Taking the spatial 
heterogeneity of R into account and neglecting terms in- 
dependent of A, the total elastic free energy is written 
as 



j dro (Tr G R : G) , 



(7) 



where uq is the number density of subchains. Inhomo- 
geneous contribution of the form R : G can be also de- 
rived from Cauchy's theory of solids bound by a central 
force We shall call Rij the quenched random 

stress although it is more directly related to quenched 
random strain in the present model. For simplicity, we 
assume that the frozen heterogeneities have a single char- 
acteristic size ^ij, which is substantially larger than the 
mesh size. After a coarse-graining on the scale ^^j, we can 
regard Rij as a spatially uncorrelated Gaussian random 
variable satisfying 



(i?.,(ro)) =0, (8) 



(9) 



The dimensionless constants (3 and p' represent the mag- 
nitudes of shear and dilatational quenched strains, re- 
spectively. 



B. Random stresses in nematic networks 

Next we consider nematic elastomers and gels. Warner 
et al. Q constructed an affine-deformation theory of 
nematic networks, by generalizing the classical theory. 
Their basic observation is that nematic chains with low 
backbone rigidity are well characterized by an anisotropic 
Gaussian conformation, elongated along the director. 
The equilibrium distribution of the end-to-end vector can 
be written in the form, 



Pe,(p) =AA'-iexp 



where a is the degree of chain anisotropy and 



1 



Qij — Qo [ ^ij ~ -^niJlj 



(10) 



(11) 



is the orientational order parameter with n being the di- 
rector. We consider a system deep in the nematic phase 
and put Qo = 1; the state of orientation is completely 
specified by the director. The coupling constant a is ex- 
pressed in terms of the parameters used in [n2l as 



^11 



(1 - + 



(12) 



Note that a does not exceed d/{d—\), the value attained 
in the anisotropic limit — > oo. An advantage of the 
affine-deformation model is that it can describe arbitrary 
crosslinking conditions; the networks can be fabricated 
either in the isotropic or the nematic phase. First we 
consider networks originally crosslinked in the isotropic 
phase, and shall treat the case of anisotropic crosslinking 
in Section V. The random stresses are now readily incor- 
porated into the original model. For the case of isotropic 
crosslinking, the initial chain conformation can be de- 
scribed by Eq. (S) , with (||) and (^) . Substituting (||) and 
( |To| ) into Eq.(^, and dropping terms independent of A 
and/or Q, we arrive at the elastic free energy. 



Fei = ^ I dr Tr 



(l + R)-A^-(l-aQ)-A-l 



(13) 



with fi — ksTuolQ/Q'). Here we subtracted a constant 
so that Fei vanishes when A = I and a — P = 0. We 
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also replaced J dro with J dr, assuming an incompress- 
ible network and imposing the local constraint, 



detA = 1. 



(14) 



We decompose the elastic free energy into proper and 
disorder contributions as Fei — + , where, by def- 
inition, the former is given by formally putting R = in 
the right hand side of ([T^), and 



Fi'^l I dvTr 



R-A^-(I-qQ)-A 



(15) 



The total free energy of the system is written as = 
Fei + Ff where Fp is the Frank free energy, for which we 
use the so-called one-constant approximation |39| , 



Fp 



K 



dr{Vn) 



(16) 



We assume that the average deformation is a uniaxial 
strain along the a;-axis, parameterized by the elongation 
ratio A, as 



A = Ae^e^ + A-i/^^^-i^l-e^e,). 



(17) 



The internal displacement is defined as the deviation 
from the average deformation, 

u^r-H-TQ. (18) 

with which the deformation tensor is expressed as 



(19) 



(here and hereafter, we imply summation over repeated 
indices i,j,k,l, and to). 

In the absence of quenched disorder, the elastic free 
energy is minimized at A = Xm and u = 0, where 



Am — 



1 + (l/d)a 



-1 {d-l)/2d 



1 - (1- l/d)c 



(20) 



is the ratio of the spontaneous elongation induced by the 
isotropic- nematic transition (see Fig.||). However, 

if the random stresses are enough strong, the long-range 
orientational order is destroyed and the ground state of 
the system is shifted to a polydomain state with A = 1, 
as we shall see. Hereafter and throughout the paper, we 
regard A as an externally controlled parameter. 



C. Effective free energy 

In this subsection, we derive the effective free energy 
in the mechanical equilibrium state under the constraint 
A = 1, and discuss its physical consequences. Substitut- 
ing Eqs.(0), & and (|l|) into Eq.(|l3|), expanding it 
with respect to Vtt and retaining a bilinear form in Vm, 
R and Q, we have 



F,j 



u=i 



dr {diUjY + 2{Rij - aQij)diUj 

-aR^jQ^j (21) 



From this we eliminate the displacement field using the 
mechanical equilibrium condition, 



SFel 

5u 



= 0, 



(22) 



and the incompressibility condition (^), which to the 
lowest order in Vu reads 



V • M = 0. 



(23) 



After a straightforward calculation following the proce- 
dure described in [0, we obtain an effective elastic free 
energy which is correct to the quadratic order in a, /?, 
and (3' , as 



F,, 



\x=i 



(I - qq) ■ [q ■ R{q) - aq ■ Q{q) 
+aRiq) : Q(-q) 



(24) 



where q — q/\q\ and Jg = /(27r) "^dq (the tilde is to put 
to express the effective nature of the free energy). The 
proper contribution to the free energy is given by 



(I - qq) ■ (q ■ Q{q)) 



(25) 



In the real space, Eq.(25) is rewritten in the form of a 
two-body long-range interaction, as 



dr / dr' 



\k{r)-d^djGi{r-r')-Qjk{r') 
+Q^Ar) ■ d,djdkdiG2{r - r') ■ Qki{r') 



(26) 



where Gn{r) {n — 1,2) are the Green functions defined 
by 



(V2)"G„(r) = ~5{r), 
G„(r ^oo) = 0. 



(27) 
(28) 



In a similar manner, the disorder part of the free energy 
is written as 



/ia dr dr' 



R^k{r) 



-A,6{r - r') -f 9,9,Gi(r - r') • Q,k{r') 



+R,,{r) ■ d,djdkdiG2{r - r') • Qki{r') 



const., (29) 
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r — r', the kernels didjGi 
^ ) decay in proportion to 



In a given direction oi R = 
and didjdkdiG2 in (^) and 

Let us discuss the physical meaning of the effective 
free energy. First we consider the disorder part, neglect- 
ing the proper elastic interaction for the moment. If we 
decompose the random stress into the dilatational part 
Rkk^ij/d and the shear (traceless) part Rij — Rkk^ijjd, 
the former makes no contribution to the free energy ( |2l| ) 
because of ( ^3|) and the tracelessness of Q. Thus, only 
the shear portion of the random stresses (whose strength 
is represented by (3) is relevant, at least in the bilin- 
ear order. It is intuitively obvious that a mere volume 
change does not create any preferential director orienta- 
tion, while anisotropic strain does. As seen from Eq.(29), 



the random stresses both locally and non-locally act on 
the director field. The classical scaling argument by Imry 
and Ma tells that arbitrary weak random stresses 
destroy the lonokrange orientation order in dimensions 
lower than four u. The domain size or the orientational 
correlation length, which we denote by ^, is determined 
by a balance between the effects of random stresses and 
Frank elasticity. The effective strength of disorder is ex- 
pressed by the dimensionless parameter. 



D = 



a- 



(31) 



A straightforward calculation reduces Eq. (^6|) to 
.2 /• r 1 



16^ 



• cos 



drj dr'- 
2(0(r)- V) + 2(0(r')-V) 



(34) 



where ip is the azimuthal angle of i? = r r' = 
|il|(cos-!/', sin-0). From the angle-dependence of the in- 
tegrand, we expect that the above free energy is mini- 
mized by a "checkered" domain configuration as depicted 
in FigJ^. Correlation in directions parallel and perpen- 
dicular to the local director is suppressed while those 
in oblique directions are enhanced. It has the following 
simple interpretation. Upon the isotropic-nematic tran- 
sition, each part of the network tends to elongate along 
the local director. The domain in the center of the figure 
pushes the top neighbor upward, pulls the left neighbor 
rightward, and so on. To reduce the mechanical con- 
flict without violating the global constraint A = 1, the 
top and left domains are reoriented perpendicular to the 
central one. This domain reconfiguration enables the I- 
N transition-induced elongation along the local director, 
despite of spatial inhomogeneity. 

The same picture holds for orientational correlation in 
three dimensions. In 3D, Eq.(E^) becomes 



According to the Imry-Ma argument, the domain size 
scales as oc Z)2/('^^4) ^^ic weak disorder regime 

D ^ 1. For a strong disorder _D > 1, we should have 
C ~ ^-R ^-iid optimization of the director field will reduce 
the disorder free energy density roughly by /ia/3. 

Next we turn to the proper elastic interaction. It 
should play only a secondary role in selecting the do- 
main size ^ because of the invariance of Eq. (^5|) against 
a change of scale q const, x q. However, it creates 
a characteristic anisotropy in the the orientational cor- 
relation. We see this first in the two dimensional case. 
In 2D, the orientational configuration is specified by the 
director's azimuthal angle 6 = 6{r), as 



n = {cos 9, sin( 



or, equivalently, 



Q 



cos 20 sin 29 
sin 26* -cos 26* 



(32) 



(33) 



F,_ 



1 



dr / dr'j^ g{n,n',R), (35) 



gin, n',R)^^-+ 4(n • n'f + (n • Rf + (n' • Rf 

o 



l^{n ■ n'){n ■ R){n' ■ R) 
15{n- Rf{n' ■ Rf, 



(36) 



where n = n(r), n' = n(r') and R — R/\R\. Correla- 
tion in the direction parallel to the director is suppressed 
as in the 2D case, which is known by observing that the 
function 



g{n, n',n) = + 2(n • n' f 



(37) 



takes its minimum when n _L n'. 

The domain reconfiguration due to the proper elastic 
interaction is suppressed by the Frank elasticity at wave- 
lengths shorter than 



^Although the original Imry-Ma argument assumes an uncorrelated random field, it is easy to see that it also holds in the 
present case, including the scaling law for the domain size. To see this, it is useful to rewrite the right hand side of Eq.(^) 
into the form /la JdrP(r) : Q(r), where the effective random field P has a long-range correlation schematically represented as 

(P,,(r)Pfe,(r')) =CflK«5(r-r')+n:,fe,|r-r'|-'], (30) 

where FT depends on the direction of r — r' but not on its magnitude. Since both fl and FT are dimensionless quantities, there 
appears no additional characteristic length that affects the Imry-Ma scaling of disorder free energy contained per domain. 
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/ K 



(38) 



Thus we have three characteristic lengthscales, ^, ^_r, and 
^c- The observed domain size ^ is typically 1 — lO^^m, 
while we estimate to be 10 nm for typical experimen- 
tal values K = 10"" J/m, fx = 10^ J/m^, and a = 1.0. 
There is a substantial gap between ^ and fc, where the 
proper elastic interaction plays a dominant role. The 
Frank free energy density fp (averaged over space) scales 
as fp ^ fei ■ < fei M"^- The domain size ^ 

can be cast into a scaling form, 

Although S is a highly non-trivial function, it can be nu- 
merically obtainable unless D is very small (or, unless 
£,/£,R is very large), as we see in Section III. We have a 
trial estimate D ^ 1 if we assume (3 ~ 0.01 and ~ 100 
nm in addition to the above values of K, fj,, and a. Of 
course, this estimate of D is quite uncertain because the 
magnitudes of /3 and should depend on the kinetics of 
the crosslinking process, quality of the solvent, etc. Our 
point here is that it is not unreasonable to have a moder- 
ately strong disorder in the presence of submicron-scale 
network heterogeneities, which is considered ubiquitous. 



D. Mechanical response 

Now we proceed to discuss the mechanical response 
during the polydomain-monodomain transition. To do 
so, it is useful to examine again the polydomain state at 
A = 1 and in 2D. The harmonic free energy ( psf ) can be 
rewritten as ISlll 



\Qiiq)\ 



Qi{q) = 2q^qyQxx{q) - (ql - ql)Qxy{q) 
= sin 2^3 Qxx{q) - cos2(p Qxy{q), 



(40) 



(41) 



where ip is the azimuthal angle of the wavevector, q = 
|<7|(cosiy9,sin93). Complementary to Qi{q) is the variable 
defined by 



h{q) = [ql - ql)Qxx{q) + 2qxqyQxyiq) 
= cos2ip Qxx{q) + sin2(^ Qxyiq)- 



(42) 



where Qa{f) (a = 1, 2) are the inverse Fourier transform 
of Qaig)- To reduce the free energy (^), there arises an 
asymmetry Qi{r)^ > Q2{i'Y ■ In the limit where /ia^ is 
much larger than the disorder and Frank contributions 
to the free energy density, we expect from (H) to have 



Qi(r)2 = -, Q2(r)2 = 0, 

which indeed is numerically confirmed p2[ 
the elastic free energy density is given by 



(45) 



In this limit, 



fei 



(46) 



as seen from (^0|). To the second order in a, it is equal to 
the free energy in the monodomain state with A = Am, 
as we can easily check by substituting (|2^) into ( [l3| ) and 
expanding it with respect to a. Thus we conclude that 
the elastic free energy change accompanied with the P-M 
transition is of 0{a^), and the macroscopic stress aver- 
aged over the region 1 < A < Am, or 



fel{X = Xrn)-.fel{X=l) 



A,, 



1 



(47) 



is a quantity of 0{a^). 

To see the origin of the soft response, it is useful to 
look at the local elastic stress tensor, which is given in 
the harmonic approximation (e|) as H] 



(48) 



Consider its variance cr?-. In the absence of random 
stresses, we have 



dr afj = /iFe; ~ J Qfj 
= t^^o? j dr{Q\-Ql) 



2j 



(49) 



which vanishes from (ffq). This means that each part 
of the system is stretched along the local director by 
1 -f a/4 4- 0(o?) ~ Am times. This local elongation, 
realized by the checkered polydomain structure, reduces 
the free energy close to its absolute minimum. 



Note that Q\(ci) and Q^isi) constitute a set of normal 
modes, and satisfy 

|Qi(g)|' + \Q2{qt = \Qxx{q)t + \Qxy{q)f. (43) 



Qi(r-)2 Q2(r)2 = Qxx{rY -f Qxy{r, 



2 _ 



4' 



(44) 



III. NUMERICAL SIMULATION 

To further study non-linear mechanical response and 
effect of random disorder, we resort to numerical simu- 
lation by the continuum model. We utilize two different 
numerical schemes, one for the polydomain state in me- 
chanical equilibrium at A = 1 and another for the P-M 
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transition and dynamical effects. A two dimensional sys- 
tem is assumed for computational advantage. All the 
simulations below are performed on a x iV square lat- 
tice with A^ — 128 unless otherwise stated. The grid spac- 
ing is chosen to be the unit of length. Periodic boundary 
conditions are imposed on n(r) and u{r), while the av- 
erage strain A is externally controlled. 



A. Polydomain state 

First we study the the polydomain state in complete 
mechanical equilibrium and with no average strain (A = 
1). To this end, we assume the harmonic free energy 
and solve the linear equations (22) and (|2^) using 
fast Fourier transform. To minimize the free energy, we 
adapted a variant of the simulated annealing method . 
The orientational order parameter is evolved according to 
a Langevin equation, 



— = r„ (I - nn) 
ot 



dF 
dn 



V 



(50) 



where r„ is a constant and 77 is a "thermal" noise satis- 
fying 



{v{r,tMr\t')}=r,l\-6{r^r')5{t~t') 



(51) 



and Gaussian statistics. The noise strength -qo is grad- 
ually reduced to zero until the end of each run. To be 
precise, we decrease 770 to zero at a constant rate in the 
former half of a run, and set 770 = in the latter half. The 
initial noise strength and the annealing rate are chosen 
so that two different initial configurations, one with ran- 
dom and another with homogeneous director field, lead 
to indistinguishable results for the macroscopic quantities 
such as correlation function, average orientation, and free 
energy densities. As a standard set of static parameters 
we choose 



(52) 



^ = 400, a = 0.2, (3 = 0.025, = 1, /sT = 4, 



for which — 0.5 and D = 0.5. We integrated Eq.(|5C 
using the Euler scheme with time increment At = 1 per 
step. A typical run consisted of 5 X 10^ time steps. Longer 
runs did not make an observable difference in the macro- 
scopic quantities. 

First we consider the orientational correlation function. 



G{R) = {Q,,{r)Q,,{r + R) 



(53) 



which is a function only of distance. We define the cor- 
relation length ^ through 



G(0) 



(54) 



For each parameter set, we took statistical average over 
20 samples. The data are shown in Fig.^. The de- 
cay of G{R) is nearly exponential for strong disorder 



and faster than exponential for weak disorder. This 
qualitative tendency agrees with previous results for the 
2D random-field XY model |||,|4[^]. The correlation 
length is a rapidly decreasing function of the effective 
disorder strength, D. The dependence is roughly expo- 
nential, also in agreement with previous results for the 
XY model In the same figure we show the de- 

pendence of ^ on /io;^ , which is the measure of elastic in- 
teraction. Although the dependence is weak, the proper 
elastic interaction has an effect of increasing the correla- 
tion length. This is related to the enhancement of corre- 
lation in directions oblique to the local director, depicted 
in Fig.^ In order to quantify the director-relative corre- 
lation, we define the function 



(55) 



H{R) = (Q,,{r)Q,j{r + \){r) ■ R) 



where U(r) is a matrix of rotation that maps n{r) to e^, 
or, explicitly, 



cos W — sm f 
sin cos t 



(56) 



By definition, H{x,0) and H{Q,y) respectively describe 
the correlation in directions parallel and perpendicular to 
the local director. The data for the standard parameter 
are plotted in Fig.|[ We see that the correlation is long- 
ranged in any specific direction, and the exponential-like 
decay in Fig.^ should be considered as a result of mu- 
tual cancellation of positive and negative correlation by 
taking the angular average. 

A real space snapshot of the order parameter field Q^y 
is also given in Fig.^. As the grayscale shows, the con- 
tour Qri-y = preferentially lies in the horizontal (x-) and 
vertical (y-) directions. This corresponds to the check- 
ered domain structure in Fig.^ (note that the grayscale 
is chosen so that the director is oblique to the horizontal 
axis in the brightest and darkest regions). More precisely, 
the checkered pattern is found on many different length- 
scales, which is a natural consequence of the fact that the 
elastic interaction energy (|2j) is scale- independent. 

An experimentally accessible way to characterize the 
anisotropic director correlation is the polarized light scat- 
tering. In a weakly inhomogeneous state, the depolarized 
(HV) light scattering intensity is given by 



Hq) = {\Q.y{q)f 



(57) 



except for a q-independent prefactor. According to 
Ref. [0 , the above formula holds even in a highly inho- 
mogeneous state, if one assumes a two-dimensional con- 
figuration (see Eq.(2) in the reference). Our numerical 
data is shown in Fig||. The intensity ( |57|) is expressed 
in terms of Qi and Q2 as I{q) = cosip {\Qi{q)\'^) + 
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sin(y3^ (|Q2(q')P), and the asymmetry Qi > Q2 explains 
the enhanced scattering on q^- and Qy- axes [plj . 

Note that the peak is located at a small but finite 
wavenumber, in contrary to what is expected from the 
non-conserved nature of the orientational order parame- 
ter. In fact, we find it to be a finite size effect, and the 
peak wavenumber shrinks to zero as the system size N 
is taken to infinity, leaving a singular minimum at the 
origin. To see this, we have computed the circularly av- 
eraged structure factor. 



Siq) 



dv ( \Q^j{Q)\ 



(58) 



for N = 64,128, and 256 systems, and found a peak in 
the region {2tt/N) < q < 2{2tt/N) in every case. The 
origin of the singular minimum at g = is explained as 
follows. Because of the periodic boundary condition on 
u, the spatial average Vm should complete vanish. This 
constraint suppresses formation of the checkered pattern 
with the check size larger than N/2. 



B. P-M transition 

Next we study the P-M transition using the non-linear 
elastic free energy (13). We found that complete mini- 
mization of the free energy takes very much computation 
time, and decided to take a more empirical approach : 
we utilize a simple dynamical model, and abandon to 
exclude non-equilibrium effects from the results. Fortu- 
nately, the stress-strain relation thus obtained is equi- 
librated to a good degree, because of fast relaxation of 
the rubber-elastic free energy. On the other hand, the 
domain structure exhibits a slow coarsening, which we 
study in the absence of external strain. 

Our dynamical model consists of a set of equations that 
describe evolution of non-conserved order parameters in 
a simplest manner, namely. 



dn dF 
- = -r„(l-nn).— , 



and 



du 
'di 



" 6u 



(59) 



(60) 



Instead of imposing the strict incompressibility condition 
( p3| ) , we penalized local volume change by adding an arti- 
ficial potential to the free energy. By taking it in the 
form F^ = ^Jdr [ao(det A - 1)^ -|- ai(det A - 1)"*] and 
choosing appropriate values of the constants ag and ai, 
we kept detA in the region [0.99,1.01] throughout the 
runs. 

We integrated (|5|) and ( |60| ) using the Euler scheme 
with r„ = 0.2, r„ = 0.02 and At = 1. A typical set of 
static parameters is same to that given by (E2). 



To prepare a polydomain state, we set site-wise ran- 
dom numbers to Q and u as the initial condition, and in- 
tegrated ( |59| ) and ( |60| ) for 5 x 10^ time steps with A = 1. 
Then we increased A at a constant rate dX/dt = Ix 10~^ 
to induce the P-M transition. To check hysteresis, fi- 
nally we decreased A back to unity at the rate dX/dt — 
-1 X 10-s. 

Plotted in Fig.^ are the scaled macroscopic elastic 
stress n~^amacro ~ fJ'^^{dfei/dX) and the mean orien- 
tation S = cos 26 = 2Qxx as functions of A. We see from 
the figure that the elastic stress is vanishingly small and 
the orientation linearly increases in the polydomain re- 
gion 1 < A < Xm{= 1.05). The stress shows a linear rise 
in the monodomain region A > Am, where the orienta- 
tion is nearly saturated to the maximum, S — I. While 
the strain-orientation curve has a small hysteresis, the 
strain-stress curve is almost completely reversible. 

The smallness of the hysteresis manifests an impor- 
tant difference between the present system and random 
anisotropy magnets under magnetic field. In the latter, 
the macroscopic orientational order is broken solely by a 
random field. In contrast, in the present model, the mon- 
odomain state is unstable to an strain-mediated director 
buckling for A < A™ , even when there is no quenched dis- 
order. This instability, which will be discussed in Section 
IV in detail, makes the P-M transition almost reversible. 

The free energy densities are also shown in Fig.||. Both 
the proper and disorder parts of the rubber-elastic free 
energy change little in the region A < A^. The latter 
curve has a slightly positive gradient. The situation is 
more subtle for the former. Its gradient is slightly pos- 
itive in the figure, and turns to slightly negative for a 
smaller disorder strength. However, in the absence of 
random stresses and at A = 1, we had four domains 
whose sizes are limited by the system size, and the do- 
main boundaries raise the elastic free energy. Because 
of this finite size effect, we cannot exactly tell the sign 
of the proper elastic stress in the macroscopic limit. We 
cannot exclude the possibility that the macroscopic stress 
completely vanishes in the limit of weak disorder. 

The strain-stress and strain-orientation curves for 
larger values of a are given in Fig.^. Each curve shows 
a sharp crossover around A — A,„(q!). The elastic stress 
in the polydomain region is vanishingly small even for 
large coupling. For any value of a studied, the changes 
of the proper and disorder elastic free energy densities, 
\fSi^ = Am) - fP{X = 1)1 and |/i?(A = A^) - fgiX = 
1)1, were smaller than 0.3 percent of fia^. We find essen- 
tially no a-dependence of the macroscopic stress in the 
polydomain region. 

Shown in Fig.^ is the histogram of the elastic free en- 
ergy contained in a lattice site. The distribution is fairly 
sharp and little changed by stretching for A < A,„, im- 
plying that the free energy is homogeneously minimized 
in the polydomain state. Real space snapshots of the do- 
main morphology is shown in Fig.^ Pinned defects are 
observed just below the threshold A = Am, while we find 
no defects remaining in the monodomain state. 
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The depolarized scattering intensity is shown in Fig.|lC|. 
It has a minimum at g = and develops four peaks at fi- 
nite wavenumbers. As we shall see in the next subsection, 
the peaks move toward the origin as the true equilibrium 
is approached and the domains coarsen. Here we con- 
centrate on the effect of stretching. The peak intensity 
first increases and then decreases as a function of A. Un- 
der stretching along the x-axis, the peaks on the Qx- axis 
are more enhanced than those on the qy- axis. The shift 
of peak wavenumber by stretching is very small and dif- 
ficult to estimate. By our choice of the stretching rate, 
the P-M transition completed in 5 x 10'^ time steps, much 
before a significant coarsening can occur. 



C. Slow structural relaxation 

Now we turn to dynamical effects. First let us discuss 
the conditions under which Eqs.(^) and ( |60|) are most 
reasonable as a model of dynamic evolution, not only as 
an artificial scheme of functional minimization. Firstly, 
Eq.(|60[) means that the velocity du/dt is proportional 
to the force —5F/6u. This applies to motion of a net- 
work in a viscous solvent p6| , where we have a straight- 
forward analogy to D'arcy's law in porous media. On 
the other hand, in dry elastomers, there arises a viscous 
stress due to intra-network friction, which is proportional 
to V{du/dt). This is not accounted for in Eq.(|6y). Thus 
we consider that the dynamic model is more appropriate 
to swollen gels than to elastomers. Secondly, Eqs.(|59|) 
and ( |60| ) neglect dynamical coupling between the order 
parameters, i.e., the non-diagonal part of the Onsager 
coefficient matrix. This does not matter if the dynamics 
of the orientational order parameter is fast and slaved to 
that of the displacement field, which we expect to be the 
case. In fact, if the constituent polymer of the gel is not 
rigid, is of the order of the viscosity of low- molecular 
weight fluids, ry. On the other hand, the friction between 
the network and solvent renders F"^ to be of the order 
of rj/P, where I is the mesh size of the network p9| , p7[ . 
Thus, the characteristic relaxation time of the strain at 
the scale of domain size is ~ (^/^)^(^ 1) times larger 
than that of Q. 

Evolution of the structure factor S{q) (as defined by 
Eq. (psf) ) is shown in Figjl^. The peak wavenumber de- 
creases and the peak intensity increases as a function of 
time. Also shown in the figure is the structure factor at 
complete mechanical equilibrium, which is obtained by 
the numerical scheme used in Section III A. The corre- 
lation length ^ and the inverse of the peak wavenumber 
go arc plotted in the middle of Figjl^. In the time region 
1 X 10"^ < i < 3 X 10^, the former is well fitted by a power 
law ^(t) cx with e — 0.23 ± 0.02, and the latter grows 
almost in parallel to the former. 

The Frank and rubber-elastic free energies are plotted 
as functions of time in Fig.|l^. While the elastic free en- 
ergy changes little after an early stage of around t = lO'^ 



time steps, the Frank energy density fp shows a slow and 
continuous decrease, which is approximately described by 
an power law fp oc t"^ in the region 1 x 10^ < i < 3 x 10^, 
with e' = 0.22 ± 0.03. 

Presently we have no explanation for the good fits of 
^(t) and fpit) by power laws. We keep ourselves to point 
out that the values of e and e' are much smaller than the 
corresponding exponents for the 2D non-conserved XY 
model without quenched disorder, which equal 0.5 and 
1.0 from a simple scaling argument |^^. The naive scal- 
ing relation e' = 2e is also broken here, which is not at all 
surprising if we consider the presence of quenched disor- 
der 1^ . We should also stress that the final equilibrium 
values of ^ and fp are finite. Preliminary study by a 
longer run without statistics finds a crossover from the 
power-law type kinetics to a slower one at t ^ 1 x 10^ 
steps. 

The above results show that the relaxation process can 
be decomposed into three characteristic stages : (i) The 
quench into the nematic phase from the isotropic phase 
produces microscopic textures, which coarsen to reduce 
both the rubber-elastic and Frank free energies. After 
the characteristic domain size reaches ^c, anisotropic do- 
main reconfiguration on this scale follows. The rubber- 
elastic free energy is almost completely minimized at 
this early stage, because of the scale-independence of 
the proper elastic interaction ([2^). (ii) The "checkered" 
domain structure further coarsens to reduce the Frank 
free energy. The domain size ^ and the peak wavelength 
27r/(3'o grow in parallel to each other, (iii) The domain 
size converges to a finite equilibrium value, while the 
anisotropic domain reconfiguration proceeds on larger 
scales (limt^oo qoii) = 0). 



IV. FLUCTUATION IN THE MONODOMAIN 
STATE 

Recall that, if there is no quenched disorder, the 
ground state of the system is the macroscopically elon- 
gated state with A — A„i and it = 0. In this state, there 
are so-called soft modes of director fluctuation, which do 
not accompany any change in the rubber elastic free en- 
ergy [p|, p^|jTl| ] . The presence of the soft modes implies 
that a homogeneous director conflguration becomes un- 
stable for A < \m\ when we compress the gel along the 
optical axis, the director "buckles" to partially cancel the 
rise of elastic free energy by compression. The result of 
the previous section means that this instability is almost 
completely soft even for large deformations. In this sec- 
tion, we look at the monodomain region A > Am and ana- 
lyze the director fluctuation modes in a harmonic level. It 
was suggested in Ref. Q that the soft fluctuations at the 
critical point A — Xm is strongly enhanced by quenched 
disorder and satisfy <^|(5n(q)p) (x q~*. Nonetheless, it 
was not fully confirmed because of a breakdown of the 
harmonic approximation at the critical point, A = A„i- 
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Also, the model used in |^ remains largely phenomeno- 
logical, and a quantitative assessment of the prediction 
is necessary. Here we extend the analysis to arbitrary 
values of A, and discuss the possibility to find disorder- 
enhanced fluctuations in practical situations. 

The director is decomposed into a homogeneous part 
and a small deviation, as 



n{r) 



S'n{r) 



(61) 



correspond to the soft modes. Similar results have been 
obtained by Olmsted fll]] for monodomain elastomers 
crosslinked in the nematic phase and without external 
strain. 

The random stresses shift the ground state to an in- 
homogeneous state, 5n = Snu and u = up. The frozen 
director deviation Snji is obtained by minimizing the to- 
tal free energy Fg/ + Fp with respect to (5n, as 



Expanding the basic free energy ( |T^ ) with respect to 5n 
and M, and then eliminating the elastic field using the 
mechanical equilibrium condition, we obtain an effective 
free energy in terms of 5n. An outline of the calculation 
is given in Appendix A. For the three dimensional case, 
the result is 



1 



KAq) + B2 qjR[M) 
^2(1 + ^1) 



A1+A2 + K'q^ 

^2(^2 + ^3) 
A1+A2 + K'q^ 



q,q.jR'j^{q) 
q^iqq: R'{q) 



, (70) 



T 



Ai(A,q) \Sn{q) 
-R[x{Q)5ni{-q) 
-Bi{X, q)qtR[M) (q ' ^n{-q 

-B2{X,q)(q-R'{q)-Sn{~q 



1 , ,2 

-A2{X,q) \q ■ Sn{q)\ 



where we have introduced a scaled Frank constant, 



(71) 



B3{X,q)[R'iq):qq][q-Sn{~q) 



A,{X,q)^X^-l- 



3a 



6-2 



3 + 4a 1 + (A^* - l)g2 ' 



A2{X,q) 
Bi{X,q) 
B2{X,q) 

Bz{X,q) 
R^,iQ) 



3a(3+4a) r^x 
3+a 3+a(2+ql) 



3a 
3-2a 



1 + (A3 - l)g2 

3a f 3+a , \3\-2 
3+^i3^ + ^ )lx 



1-1- (A3 
X^qx 



)ql 



1 + (A3 

( 3+a 



3-2q 



_l + (A3-l)gr 

AikAjiRki{q), 



(63) 

(64) 

(65) 

(66) 

(67) 
(68) 



At the critical point A — Am, the quantity Ai + A2 
appearing in ( |70| ) vanishes for q _L e^;. Hence, in the 
long wavelength limit (7 — > 0, we have Snfi(q) oc q"^ in 
the soft directions q ± and q \\ e^- This means a 
divergence of the real space amplitude (|nfl(r)p) and 
(62)b reakdown of the harmonic approximation, as pointed 
out in Ref. Q. Severer is the divergence of the frozen 
elastic field uji, which is related to Snu through the 
mechanical equilibrium condition (A2). At A = Am, it 
behaves as (|Mfl(q)p) cx q~^ in the soft directions, and 
(|Mi?(r)p) diverges. However, these divergences disap- 
pear for A > Am, where the excess stretching acts as 
a stabilizing field. Now we consider this region. The 
condition for the harmonic approximation to be valid is 
{\5nji{r)\'^) < 1, which implies (|Vufi(r)p) < a^. To 
assess the condition by order estimate, we concentrate on 
the plane q _L e^;, from where arises the most significant 
contribution to the real space amplitude. 



\SnR{r)\' 



(system's volume) Jg 



\6nniq)n. (72) 



which is correct to the bilinear order in 5n and R (we 
neglect terms independent of Sn). 

In the absence of quenched disorder, the integrand 
in (^ is of the form ^{Ai\ + A2qq) : 5n(q)5n{-q). 
It is convenient to introduce two unit vectors ei — 
q X ex/\q x ex\ and 62 = ei x ex |^, with which the 
integrand becomes 



On that plane, the strongest q-dependence of Snji{q) 
comes from a factor {A + K'q"^)'^ , where 



A = A{X)=[A,{X,q)+A2{X,q)]\ 



(73) 



lA,\e, 



■ Sn{q)\ 



^(^1+^2(1 



|e2-^n(q)| . (69) 



is the measure of excess stretching. Note that A cx A— Am 
for A — Am <C 1- A saddle-point approximation around 
the plane yields 

,^2f3„2 



At A = Am, the coefficient Ai takes its minimum value 
for q II Ex, while Ai+A2{l—q^) is minimized and vanishes 
both on the line q || e^, and in the plane q J- e^- These 



C3 



(A 
K'q 



2 

■max 



;^/3/2 



K'q, 



2 

max 



(A « K'ql 
(A » K'q 



3/2 



(74) 



^ ) 

max I 
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where qmax is the upper cutoff wavenumber. We may 
roughly identify ^T^jqmax with the network mesh size, 
I ~ {ksT / nY^^ . Using typical experimental values 
H = 10^ J/m3, K = 10-" J/m, a = 1.0 and T = 300 K, 
we have = 10 nm and K'q^^^ ~ 1. For A < K'q^^^ 
and except in the close vicinity of the criticality, A = 0, 
the amplitude only weakly depends on the elongation ra- 
tio. In this region, the harmonic approximation is valid 
if and only if 



« 1, 



(75) 



10^ nm and (3 - 0.01 



which is satisfied when we put 
as a trial. 

The director exhibits a thermal fluctuation around the 
inhomogeneous ground state. Its amplitude is not af- 
fected by the quenched randomness, at least within the 
harmonic calculation. The total fluctuation amplitude is 
given by 



knTX 



1 



Pl^\q) = {\e,-6ni,{qf), 



(76) 
(77) 
(78) 



where a — 1,2, and P^^ and P^"'' arc the thermal 
and frozen contributions, respectively. Let us com- 
pare the two contributions. To be explicit, we com- 



(a) 



pare p!j?\q) and P^'{q) on the plane q _L e^. There 



>(2) 



(2) (2) 

the ratio / Pj, is controlled by a factor of the form 
Ac/(A + K'q^), where 



A, = 



(79) 



defines a crossover point. In the region A < Ac, the dis- 
order part of the fluctuation dominates the thermal part 
at long wavelengths. We estimate Ac ~ 10"^ using the 
above mentioned values, for which the crossover length 
^jK^jK^c. ^ 10^ — 10^ nm is around or below the wave- 
length of visible light. The amplitudes for a = 0.2 and 
A/A„i = 1.001 are plotted in Fig.|l^. We see that the 
anisotropy of the scattering pattern is not much affected 
by the quenched disorder. 

The director fluctuation amplitudes are closely related 
to polarized light scattering intensity By compar- 

ing experimental results to the above calculation, we may 
extract information on the network heterogeneity. In par- 
ticular, if the macroscopic orientation is not saturated in 
the monodomain state, it means the presence of large- 
scale quenched strains that do not meet the condition 
([zsl). For instance, in an optical study of a swollen mon- 
odomain gel by Chang et al. |5^], speckles on the few- 
micrometer scale are observed, which is attributed to 
heterogeneities as we consider here. We hope that the 
scattering intensity will be measured as a function of ap- 
plied strain. Another origin of large scale heterogeneity 
will be discussed in the next section. 



V. EFFECT OF CROSSLINKING CONDITION 



In this section, we consider the case of anisotropic 
crosslinking. Melts of nematic polymers often exhibit 
long-lived polydomain textures after a quench from the 
isotropic phase |]5l| , ^ , ^ . The size of the domains is 
macroscopic and typically of micron order. When such 
a melt is crosslinked, its non-uniform orientation is im- 
printed into the network. We denote the initial configu- 
ration by Qo('")- The extended affine-deformation theory 
prescribes the elastic free energy, 



Fe; = I / dr Tr 



(l + aoQo)-A^-(l-aQ)-A-l 



(80) 



where oiq is expressed in terms of the parameters used 
in 1 12 1 as 



(i/d)^ii + (1 - ijdYx 



n) 



or equivalently, (Xq — a/[l — (1 — 2/(i)a]. Note that the 
above free energy is obtained by just formally replacing 
R with aoQo in the free energy ( |l3[) . Thus, the initial 
texture field Qg provides a source of quenched disorder. 
An effective disorder strength that corresponds to ( |3l| ) 
can be defined by 



(82) 



where is the correlation length of the initial texture. 
If we set ^0 = 1/im and a — 0.1, we have a very large 
number D ^ i^o/^c)'^ — 10^ ^ 10^. Since the orien- 
tational order is predominantly affected by this strong 
disorder, we do not take into account other mesoscopic 
sources of quenched stresses, which is legitimate as a first 
approximation. 

We have simulated the P-M transition in the following 
way. To generate the initial configuration, we mimicked 
the phase ordering kinetics of nematic polymer solutions 
by numerically solving the equation. 



— --F (I 

dt~ 



nn 



dFp 
dn 



(83) 



Taking a sitewise-random director configuration as the 
initial condition, we integrated Eq.(|3|) over 50 time steps 
(with F„ = 0.2 and At = 1) to generate Qo(''")- After 
crosslinking the system by adding F^i to the free energy 
and setting A = 1 and tt = 0, we integrated Eqs.([59|) and 
( |60| ) for 1 X 10'* time steps to equilibrate the system. The 
mechanical response was studied in just the same way as 
described in Section III. 

Fig.|l^ shows the strain-stress and strain-orientation 
curves for a = 0.2,0.4,0.8, and 1.2. The strain-stress 
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curve bends at a value of A where the director is yet far 
from ahgned. For a > 0.8, the gradient of the strain- 
stress curve has a non-monotonic dependence on the 
strain, and is smallest at an intermediate value of A. The 
strain-orientation curve shows only a gradual crossover to 
the monodomain state, especially for larger values of a. 
The slope of the scaled elastic stress fJ.~^(Jmacro is roughly 
independent of a in the vicinity of the point A = 1. 

Evolution of the director texture during the P-M tran- 
sition is shown in Fig, 14, and the distribution of the 
elastic free energy in Fig.^. At A = 1 the director 
texture is almost same as that just before crosslinking, 
or Q(r) = Qo{r). This is reflected in the extremely 
homogeneous free energy distribution. External strain 
strongly dehomogenize the distribution, and the peak is 
continuously broadened as we increase A toward the mon- 
odomain region. 



VI. DISCUSSION AND SUMMARY 

We have studied polydomain nematic networks from 
two aspects, namely, (i) breaking of long-range orienta- 
tional order by frozen internal stresses and (ii) a non-local 
inter-domain interaction arising from strain-orientation 
coupling. The mechanical response is controlled by the 
latter if the quenched disorder is of moderate strength 
(or, if D ^ 1). In this case, the proper elastic interac- 
tion reorganizes the polydomain structure so that local 
elongation along the director is achieved everywhere in 
the system. The resulting structure contains the check- 
ered correlation pattern on various scales, which produces 
the "four-leaf clover" pattern in the depolarized scatter- 
ing intensity. Upon stretching, the director and the lo- 
cal strain axis coincidently rotate toward the direction 
of macroscopic extension, and thus the elastic free en- 
ergy keeps almost constant until a complete alignment is 
attained. The change in elastic free energy accompany- 
ing the P-M transition is analytically estimated to be of 
0{a^). This result does not depend on a specific model of 
non-linear elasticity. In fact, we obtained it by harmonic 
expansion of the elastic free energy, which is unique from 
symmetry Q . Numerical simulation reveals a more com- 
plete softness, and we find essentially no a-dependence of 
the average macroscopic stress (p7|). We cannot exclude 
the theoretical possibility that the P-M transition is ex- 
actly soft in the weak disorder limit. We may say that 
there are mechanical quasi- Goldstone modes, which are 
distinguished from the genuine Goldstone modes of fluid 
nematic liquid crystals in that there is an anisotropic cor- 
relation even in the absence of external field and that the 
quenched disorder selects a characteristic lengthscale. 

We have discussed two sources of quenched disorder. 
One is the random stress due to residual heterogeneous 
strains at the moment of crosslinking, considered ubiqui- 
tous in rubbery networks. The anisotropic (shear) part 
of random stresses act on the orientational order, both 



locally and non-locally. The macroscopic domain size ob- 
served in experiments can be explained if there are frozen 
heterogeneities of a reasonably small magnitude (e.g. 1 
percent in strain) and a size somewhat larger than that 
of individual network meshes (e.g. 10^ nm). 

A different viewpoint is taken in previous theo- 
ries pO| , pl 22 1, where a random molecular field operating 
at crosslinks is assumed to be the source of disorder. The 
random field hypothesized there has a small correlation 
length roughly equal to the distance I between crosslinks. 
The ratio is a large number (~ 10^), which means a 
weak effective disorder. 

Currently we know of no firm experimental indica- 
tion of the disorder strength. A possible method of its 
estimate is to observe director fiuctuation in the mon- 
odomain state. We have calculated thermal and disorder 
contributions to the fiuctuation amplitude, which is pro- 
portional to the polarized light scattering intensity. The 
intensity diminishes as we stretch the network, and there 
is a region of macroscopic strain A where the disorder con- 
tribution dominates over the thermal one. The width of 
the region and the absolute value of the intensity should 
inform us the order of the disorder strength. The thermal 
and quenched contributions could be separately analyzed 
by use of dynamical light scattering (DLS). Indeed, DLS 
has been successfully used to decompose the two kinds of 
density fiuctuation in gels ||^ . It is hoped that a similar 
method will be developed for orientation fluctuation in 
the present system. 

By crosslinking the network in the nematic phase and 
in the course of phase ordering, we obtain another kind of 
quenched stresses. The polydomain texture of the liquid- 
crystalline polymer melt is almost completely frozen by 
crosslinking, if its characteristic size is larger than ^c- 
The memory of the initial macroscopic texture makes the 
mechanical response non-soft. Spatial distribution of the 
elastic free energy is strongly dehomogenized by applied 
strain, in contrast to the case of isotropic crosslinking. 

The influence of crosslinking conditions has little been 
discussed in previous studies of the P-M transition, ex- 
cept for a few experimental papers [ p^Jl7| . Kiipfer and 
Finkelmann |l^] studied both isotropic and anisotropic 
crosslinkings under external stress of various magnitudes. 
Fig. 8 in the reference shows that polydomain networks 
crosslinked in the nematic phase are harder than those 
prepared in the isotropic phase. Another example of soft 
and non-soft P-M transitions is given in Ref . [jl^ , where it 
is stated that some of their samples were prepared above 
the isotropic- nematic transition temperature of the melt, 
while the others are crosslinked below it. Unfortunately, 
they do not explicitly state the crosslinking condition for 
each stress-strain curve. We wish further effort in this di- 
rection to be made in the future, especially to flnd more 
evidence of vanishing macroscopic stress. 

A remark should be made in relation to this. We have 
assumed that the quenched heterogeneities have meso- 
scopic sizes in the case of isotropic crosslinking. However, 
if the network is crosslinked in poor solvents or near the 
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spinodal line, the heterogeneities can be macroscopic and 
cause strong effective disorder. Therefore, the mechani- 
cal response should be discussed in terms of the size of 
the heterogeneity, not only on the phase where the gel is 
fabricated. Another problem in interpretation of strain- 
stress data arises from slowness of dynamical relaxation. 
A recent dynamic measurement by Clarke and Terent- 
jev Q strongly suggests that the stress level will be sub- 
stantially lowered in the final equilibrium state, which is 
not reachable on a practical timescale. It might be pos- 
sible that a soft equilibrium P-M transition is masked 
behind a stress plateau of a sizable height, which is re- 
ported in earlier studies p3| , p^ . 

We have studied dynamical relaxation after a quench 
from the isotropic phase. The structure factor develops 
a peak at a finite wavenumber, which goes to zero as the 
true equilibrium is approached. Both the inverse peak 
wavenumber and the correlation length show a power- 
law type growth in an intermediate stage, while the elas- 
tic free energy is almost completely minimized in an early 
regime of the coarsening process. 

Some of the experimentally observed features of the 
"four-leaf clover" scattering pattern have been repro- 
duced in the present work. Firstly, we propose that 
the finiteness of the observed peak wavenumber is ex- 
plained by the slow relaxation. The experimental peak 
wavenumber does not change during the P-M transi- 
tion |16|. Together with our simulation result, it sug- 
gests that the coarsening is very slow and does not occur 
in the timescale of observation. Further experimental 
study of structural relaxation in conjunction with stress 
relaxation would be informative to check this point. Sec- 
ondly, the peak intensity increases and then decreases as 
we stretch the gel. Qualitatively the same monotonic be- 
havior is reported in the experiment. The initial increase 
is due to a sharpening the peak, which is partially under- 
stood by the fact that the director fluctuation at A = A™ 
is soft only on a plane and a line in the q-space. 

We close by listing some open questions, (i) We did 
not answer whether the long-range order is destroyed by 
an arbitrarily weak disorder under no external stress or, 
equivalently, when the average strain A is not externally 
constrained. A shift of the ground state from the mon- 
odomain (A = Am) to polydomain (A = 1) states should 
occur, either gradually or abruptly, as we increase the dis- 
order strength from zero. Probably this problem is not 
of practical importance because of a small but finite hys- 
teresis and the slow dynamics, (ii) Stretching-induced- 
anisotropy of the depolarized scattering pattern as we 
numerically find is contrary to the experimental observa- 
tion. We may suggest an effect of spatial dimensional- 
ity. In three dimensions there are three Frank constants, 
whose relative strengths may affect the anisotropy. Ex- 
perimental investigation of 3D domain structure would 
be informative, (iii) Much remains to be done for un- 
derstanding dynamical relaxation to the final equilib- 
rium state. In theoretical part, the origin of the appar- 
ent power law is yet unknown. Dynamic equations for 



dry elastomers are to be constructed, taking the intra- 
network friction account. In numerical part, late stages 
of the relaxation process is left unexplored. Stress relax- 
ation for strong quenched disorder and after stretching 
should be addressed to make a comparison to experi- 
ment. As these necessitate extensive computation, we 
leave them for future work. 
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APPENDIX A: EFFECTIVE FREE ENERGY IN 
THE MONODOMAIN STATE 



Here we sketch the derivation of Eq.(62). Substituting 
Eqs.(0), dH) and (|l9|) into Eq.(|l|), we have 



dr 



+C.ijLki{diUk){djUi) + K{diUif 



(Al) 



where Cy = (4; + Rki)Mk^ji and Ly = S,j - aQtj = 
{l+a/3)Sij—aninj . In the third term of the integrand we 
have replaced Cy and with their spatial averages as 
the deviations will contribute only to higher order terms 
in the effective free energy. The last term is added by 
hand to temporarily relax the incompressibility condi- 
tion (|2^), which is recovered by taking the limit k oo 
afterwards. The condition of mechanical equilibrium ( [2^ ) 
can be written as 

di{CikLji) + CijLkididjUk + nd^djUj = (A2) 

Taking the incompressible limit k — s- cx), we have 

1 



u{q) = 



C : qq 



L 1 : qq 

where g is an auxiliary variable defined by 
5i(r) = dj{CjkLik). 



q-l-'-g{q) 



(A4) 



Substituting (A3) into ( |Al| ), we obtain an effective free 
energy. 



Fei = ^ I drC : L 



C : qq 



L ^-.qq 



-flr(q) • L -gi-q) 



(A5) 
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We arrive at Eq.(|6^) by putting (|T^) into dj, ( |6ll) i nto 
Lij, and the resulting expressions into (|A4|) and (|A5D . 



[1] K. Dusek and W. Prins, Adv. Polym. Sci. 6, 1 (1969). 
[2] For a recent review, see: M. Shibayama, Macromol. Chem. 

Phys. 199, 1 (1998). 
[3] E. Mendes, P. Lindner, M. Buzier, F. Bone and J. Bastide, 

Phys. Rev. Lett. 66, 1595 (1991). 
[4] J. Bastide, L. Leibler and J. Prost, Macromolecules 23, 

1821 (1990). 
[5] A. Onuki, J. Phys. II (France) 2, 45 (1992). 
[6] For a review, see: J. Bastide and S. J. Candau, in Physical 

Propertires of Polmeric Gels, edited by J. P. Cohen Addad, 

(Wiley New York, 1996), Chapter 5. 
[7] P. G. de Gennes, C. R. Acad. Sci. B281, 101 (1975). 
[8] M. Warner, K. P. Gelling, and T. A. Vilgis, J. Chem. Phys. 

88, 4008 (1988). 
[9] L. Golubovic and T. C. Lubensky, Phys. Rev. Lett. 63, 

1082 (1989). 

[10] M. Warner, P. Bladon, and E. M. Terentjev J. Phys. II 

France 4, 93 (1994). 
[11] P. D. Olmsted, J. Phys. II France 4, 2215 (1994). 
[12] For a review, see : M. Warner and E. M. Terentjev, Prog. 

Polym. Sci. 21, 853 (1996), and references therein. 
[13] J. Schatzle, W. Kaufhold, and H. Finkelmann, Macromol. 

Chem. 190, 3269 (1989). 
[14] J. Kiipfer and H. Finkelmann, Macromol. Chem. Phys. 

195, 1353 (1994). 
[15] G. H. F. Bergmann, H. Finkelmann, V. Percec, and M. 

Zhao, Macromol. Rapid. Commun. 18, 353 (1997). 
[16] S. M. Clarke, E. M. Terentjev, I. Kundler, and H. Finkel- 
mann, Macromolecules 31, 4862 (1998). 
[17] E. R. Zubarev, R. V. Talroze, T. I. Yuranova, N. A. Plate, 

and H. Finkelmann, Macromolecules 31, 3566 (1998). They 

obtained both soft and non-soft P-M transitions, as seen 

in their Figs. 7 and 6, respectively. 
[18] S. M. Clarke and E. M. Terentjev, Phys. Rev. Lett. 81, 

4436 (1998). 

[19] Early experiments on polydomain elastomers are reviewed 
in : H. R. Brand and H. Finkelmann, in Handbook of Liq- 
uid Crystals Vol. 3, edited by D. Demus et al. (Wiley- VCH, 
New York, 1998), Chapter V. 

[20] A. ten Bosch and L. Varichon, Macromol. Theor. Simul. 3, 
533 (1994). 

[21] S. M. Clarke, E. Nishikawa, H. Finkelmann, and E. M. 

Terentjev, Macromol. Chem. Phys. 198, 3485 (1997). 
[22] S. V. Fridrikh and E. M. Terentjev, Phys. Rev. Lett. 79, 

4661 (1997); Phys. Rev. E. 60, 1847 (1999). 
[23] L. D. Landau and E. M. Lifshitz, Theory of Elasticity, 3rd 

Edition (Pergamon, Oxford, 1986). 
[24] K. H. Lau and W. Kohn, Surf. Sci. 65, 607 (1977). 
[25] V. I. Marchenko and A. Y. Parshin, Sov. Phys. JETP 52, 

129 (1981). 

[26] J. W. Cahn, Acta MetaU. 9, 795 (1961). 

[27] A. Onuki, J. Phys. Soc. Japan. 58, 3065 (1989); ibid. 3069 



(1989). 

[28] T. Tanaka, S. -T. Sun, Y. Hirokawa, S. Katayama, J. 

Kucera, Y. Hirose, and T. Amiya, Nature 325, 796 (1987). 
[29] For a review, see: A. Onuki, Adv. Polym. Sci. 109, 63, and 

references therein. 
[30] M. Goulian, R. Bruinsma, and P. Pincus, Europhys. Lett. 

22, 145 (1993). 

[31] N. Uchida and A. Onuki, Europhys. Lett. 45, 341 (1999). 

[32] N. Uchida, Phys. Rev. E 60, R13 (1999). 

[33] P. G. de Gennes, Scaling Concepts in Polymer Physics, 
(Cornell University Press, Ithaca, 1980). 

[34] S. Alexander, J. Phys. (France) 45, 1939 (1984). 

[35] For a work in this direction, see : S. Panyukov and Y. 
Rabin, Macromolecules 29, 7960 (1996). 

[36] P. J. Flory, Principles of Polymer Chemistry (Cornell Uni- 
versity, Ithaca, 1953). 

[37] A. E. H. Love, A Treatise on the Mathematical Theory of 
Elasticity (fourth edition), (Cambridge University Press, 
Cambridge, 1927). 

[38] P. J. Flory, Trans. Faraday Soc. 56, 722 (1960). This pa- 
per derives the free energy of a network crosslinked in 
two stages. If we assume a two-stage crosslinking with null 
crosslinks introduced in the first stage and succeeding de- 
formation R before the second one, we obtain Eq.(P). 

[39] P. G. de Gennes and J. Prost, The Physics of Liquid Crys- 
tals, 2nd Edition (Oxford University Press, Oxford, 1993). 

[40] Y. Imry and S. -K. Ma, Phys. Rev. Lett. 35,1399 (1975V 

[41] The elastic stress tensor for the full non-linear model (113) 
is given by 

aij = ^{Sik + Rik)i5lm - aQlm)^jlJ^km- (A6) 

Expanding it with respect to Vu and subtracting a hydro- 
static term proportional to 5ij, we obtain Eq.([48|). 

[42] W. H. Press, A. A. Teukolsky, W. T. Vetterling, and B. 
P. Flannery, Numerical Recipes in C, 2nd edition, (Cam- 
bridge University Press, Cambridge, 1992). 

[43] B. Dieny and B. Barbara, Phys. Rev. B 41, 11549 (1990). 

[44] M. J. P. Gingras and D. A. Huse, Phys. Rev. B 53, 15193 
(1996). 

[45] Y. -K. Yu, P. L. Taylor, and E. M. Terentjev, Phys. Rev. 
Lett. 81, 128 (1998). 

[46] To be strict, flow of the solvent should be taken into ac- 
count j2^. The relation du/dt oc 5F/5u holds in the limit 
where the solvent viscosity is infinitely large. A similar sim- 
plified dynamic equation has been used, for instance, for 
describing swelling dynamics in : T. Hwa and M. Karder, 
Phys. Rev. Lett. 61, 106 (1988). 

[47] J. L. Harden, H. Pleiner, and P. A. Pincus, J. Chem. Phys. 
94, 5208 (1990). 

[48] A. J. Bray Adv. Phys. 43, 357 (1994). 

[49] Even in the absence of quenched disorder, the simple scal- 
ing result e' = 2e is ap par ently broken due to a strong finite 
size effect in 2D. See [^ and M. Zapotocky, P. M. Gold- 
bart, and N. Goldenfeld, Phys. Rev. E 51, 1216 (1995). 

[50] C. -C. Chang, L. -C. Chien, and R. B. Meyer, Phys. Rev. 
E 56, 595 (1997). 

[51] T. Takebe, T. Hashimoto, B. Ernst, P. Navard, and R. S. 
Stein, J. Chem. Phys. 92, 1386 (1990). 

[52] F. Elias, S. M. Clarke, R. Peck and E. M. Terentjev, Eu- 
rophys. Lett. 47, 442 (1999). 



14 



[53] M. Doi and S. F. Edwards, The Theory of Polymer Dy- 
namics, (Oxford University Press, Oxford, 1986), Section 



10.6 and references therein. 




FIG. 1. Schematic illustration of disordered network 
structure. When the nematic order is introduced, the di- 
rector is preferentially oriented along the extensional axes of 
frozen network strains. 



isotropic nematic 





FIG. 3. Preferred local director configuration. Ellipses in- 
dicate the anisotropy of local strain. The "checkered" struc- 
ture enables each domain to elongate along the local director, 
without violating the global constraint A = 1. Note that the 
proper elastic interaction does not select the characteristic 
size of the pattern. 




FIG. 4. Left: Correlation function G{R) for different disorder strengths. Middle: Correlation length versus disorder strength. 
In the left and middle plots we fix ij,a^ = 16. Right: Correlation length versus strength of elastic interaction. The disorder 
strength is fixed {D = 0.5). 
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FIG. 5. Loft: Director-relative correlation function H{R). Plotted in the region ^ < -R < 5^. Middle; Depolarized light 
scattering intensity {\Qxyi(i)f'j . Right: Snapshots of the orientational order parameter field Qxy{r) = sin 20. The value of Qxy 
is positive in white regions and negative in black regions. 
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FIG. 6. Top: Strain-stress (A vs. f^~^(Tmacro) and 
strain-orientation (A vs. S) curves for a = 0.2 (A^ = 1.051). 

Diamonds (<*) and crosses (+) are for the loading and unload- 
ing processes, respectively. Bottom: Free energy densities 
(for the loading process). 
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FIG. 7. Strain-stress and strain-orientation curves for 
a = 0.2, 0.4, 0.8 and 1.2 from the left to right. Correspond- 
ing values of Am are 1.05, 1.11, 1.24 and 1.41, respectively. 
The parameters iia^ = 16 and D = 0.5 are common to all 
cases. 
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FIG. 8. Histogram of the elastic free energy contained in 
a site. Distributions for cases A = 1,1.03,1.05 and 1.07 are 
shown. The first three are indistinguishable from each other. 
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FIG. 9. Real space snapshot of the field Qxyi'p)- 
I(q) 




FIG. 10. Depolarized scattering intensity at A = 1,1.03, and 1.05 from the left to right. Statistical average over 20 samples 
is taken for each case. Note the difference of intensity scales. 



17 











- 






0.22 




j: 


^ 








iDDOffl) imm Mmm 





1000 10000 100000 

tmie(steps) 



10000 100000 
time (steps) 



FIG. 11. Loft: Circularly averaged structure factor S{q) at t = 2 x 10^,8 x 10*, 32 x 10* steps and in mechanical equilibrium. 
Statistical average over 20 samples is taken. Middle: Evolution of the correlation length ^ and inverse peak wavenumber q^^. 
Right: Dynamic relaxation of Prank and elastic free energies. Plotted are data from 20 individual runs. 
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PIG. 12. Square amplitudes of director fluctuation (in ar- 
bitrary unit). P^^ and P^-' are multiplied by a common 
prefactor. Shown is the region -30 < K'^^'^qx < 30 and 
-30 < /sT'^/^q • 62 < 30. 



FIG. 13. Strain-stress and strain-orientation curves for 
the case of anisotropic crosslinking. Plotted for cases 
a = 0.2, 0.4, 0.8, and 1.2 from the left to right, with /ua^ = 16 
in common. 
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FIG. 14. Snapshots of the field Qxyir) for the case of 
anisotropic crosslinking. The director texture of the nematic 
liquid just before crosslinking is retained at A = 1. 
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